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Cell motion inside dense tissues governs many biological processes, including embryonic development and 
cancer metastasis, and recent experiments suggest that these tissues exhibit collective glassy behavior. To make 
quantitative predictions about glass transitions in tissues, we study a self-propelled Voronoi (SPY) model that 
simultaneously captures polarized cell motility and multi-body cell-cell interactions in a confluent tissue, where 
there are no gaps between cells. We demonstrate that the model exhibits a jamming transition from a solid-like 
state to a fluid-like state that is controlled by three parameters: the single-cell motile speed, the persistence 
time of single-cell tracks, and a target shape index that characterizes the competition between cell-cell adhesion 
and cortical tension. In contrast to traditional particulate glasses, we are able to identify an experimentally 
accessible structural order parameter that specifies the entire jamming surface as a function of model parameters. 
We demonstrate that a continuum Soft Glassy Rheology model precisely captures this transition in the limit of 
small persistence times, and explain how it fails in the limit of large persistence times. These results provide 
a framework for understanding the collective solid-to-liquid transitions that have been observed in embryonic 
development and cancer progression, which may be associated with Epithelial-to-Mesenchymal transition in 
these tissues. 


Recent experiments have revealed that cells in dense bio¬ 
logical tissues exhibit many of the signatures of glassy ma¬ 
terials, including caging, dynamical heterogeneities and vis¬ 
coelastic behavior ||2] [31 [33 [43] H?). These dense tissues, 
where cells are touching one another with minimal spaces in 
between, are found in diverse biological processes including 
wound healing, embryonic development, and cancer metasta¬ 
sis. 

In many of these processes, tissues undergo an Epithelial- 
to-Mesenchymal Transition (EMT), where cells in a solid¬ 
like, well-ordered epithelial layer transition to a mesenchy¬ 
mal, migratory phenotype with less well-ordered cell-cell in¬ 
teractions Eiiiia, or an inverse process, the Mesenchymal- 
to-Epithelial Transition (MET). Over many decades, detailed 
cell biology research has uncovered many of the signaling 
pathways involved in these transitions iniiMi, which are im¬ 
portant in developing treatments for cancer and congenital dis¬ 
ease. 

Most previous work on EMT/MET has focused, however, 
on properties and expression levels in single cells or pairs of 
cells, leaving open the interesting question of whether there is 
a collective aspect to these transitions: Are some features of 
EMT/MET generated by large numbers of interacting cells? 
Although there is no definitive answer to this question, sev¬ 
eral recent works have suggested that EMT might coincide 
with a collective solid-to-liquid jamming transition in biolog¬ 
ical tissues miiMiiMiiiii. Therefore, our goal is to develop 
a framework for jamming and glass transitions in a minimal 
model that accounts for both cell shapes and cell motility, in 
order to make predictions that can quantitatively test this con¬ 
jecture. 

Jamming occurs in non-biological particulate systems (such 
as granular materials, polymers, colloidal suspensions, and 
foams) when their packing density is increased above some 
critical threshold, and glass transitions occur when the fluid is 
cooled below a critical temperature. Over the past 20 years 


these phenomena have been unified by “jamming phase dia¬ 
grams” 12911551 . 

Building on these successes, researchers have recently used 
self-propelled particle (SPP) models to describe dense biolog¬ 
ical tissues ll4l|T5l|20l|45l|48l|5ll. These models are similar 
to those for inert particulate matter - cells are represented as 
disks or spheres that interact with an isotropic soft repulsive 
potential - but unlike Brownian particles in a thermal bath, 
self-propelled particles exhibit persistent random walks. 

SPP models typically exhibit a glass transition from a dif¬ 
fusive fluid state to an arrested sub-diffusive solid that is con¬ 
trolled by (1) the strength of self-propulsion |[T5ll20ll35l and 
(2) the packing density (|) ||6l[T3][14]|20||35l . Just like in ther¬ 
mal systems, a jamming transition occurs at a critical packing 
density (|)g, but this critical density is altered by the persistence 
time of the random walks il[l3l[Il[20l[35l. 

During many biological processes, however, a tissue re¬ 
mains at confluence (packing fraction equal to unity) while 
it changes from a liquid-like to a solid-like state or vice-versa. 
Eor example, in would healing, cells collectively organize to 
form a ‘moving sheet’ without any change in their packing 
density ll26l . and during vertebrate embryogenesis mesendo- 
derm tissues are more fluid-like than ectoderm tissues, despite 
both having packing fraction equal to unity im. 

Recently, Bi and coworkers |[7l have demonstrated that the 
well-studied vertex model for 2-D confluent tissues ||8]n2]l22] 
exhibits a rigidity transition in the limit of zero 
cell motility. Specifically, the rigidity of the tissue vanishes 
at a critical balance between cortical tension and cell-cell ad¬ 
hesion. An important insight is that this transition depends 
sensitively on cell shapes, which are well-defined in the ver¬ 
tex model. While promising, vertex models are difficult to 
compare to some aspects of experiments because they do not 
incorporate cell motility. 

In this work, we bridge the gap between the confluent tis¬ 
sue mechanics and cell motility by studying a hybrid between 
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FIG. 1. Analysis of glassy behavior. (A) The mean-squared displacement of cell centers for = \ and vq = 0.1 and various values of 
Pq (bottom to top: po = 3.5,3.65,3.7,3.75,3.8,3.85) show the onset of dynamical arrest as po is changed indicating a glass transition. The 
dashed lines indicate a slope of 2(ballistic) and 1 (diffusive) on log-log plot. (B) The self-intermediate scattering function at the same values of 
po shown in (A) shows the emergence of caging behavior at the glass transition. (C) The effective self-diffusivity as function of po at vq = 0.1. 
At the glass transition D^ff becomes nonzero. (D)The cell displacement map in SPV model for a fluid state very close to the glass transition 
(po = 3.8, Vq = 0.1 and D,- = 1) over a time window t = 10^ corresponding to the structural relaxation at which Fi(r) ~ 1/2. 


the vertex model and the SPP model, that we name Self- 
Propelled-Voronoi (SPV) model. A similar model was in¬ 
troduced by Li and Sun 1^ . and cellular Potts models also 
bridge this gap ClEol, although glass transitions have not 
been carefully studied in any of these hybrid systems. 

I. THE SPV MODEL 

While the vertex model describes a confluent tissue as a 
polygonal tiling of space where the degrees of freedom are 
the vertices of the polygons, the SPV model identihes each 
cell only using the center (r,) of Voronoi cells in a Voronoi 
tessellation of space (Dirichlet domains) ||27l. The observa¬ 
tion that Voronoi tessellations can describe cellular patterns 
in epithelial tissues was hrst proposed by Honda llTD . For 
a tissue containing V-cells, the inter-cellular interactions are 
captured by a total energy which is the same as that in the 
vertex model. Since the tessellation is completely determined 
by the {r;}, the total tissue mechanical energy can be fully 
expressed as E = Fdr;}): 

N 

E = Y,[KA{A{ri)-Aof+Kp{P{ri)-PQf] . ( 1 ) 

1=1 

The term quadratic in cell area A(r,) results from a combi¬ 
nation of cell volume incompressibility and the monolayer’s 
resistance to height fluctuations ll22l . The term involving 
the cell perimeter P(r,) originates from active contractility 


of the acto-myosin sub-cellular cortex (quadratic in perime¬ 
ter) and effective cell membrane tension due to cell-cell ad¬ 
hesion and cortical tension (both linear in perimeter). This 
gives rise to an effective target shape index that is dimension¬ 
less: po = Po/^/Aq. Ka and Kp are the area and perimeter 
moduli, respectively. For the remainder of this manuscript we 
assume po is homogenous across a tissue, although heteroge¬ 
neous properties are also interesting to consider ||5^ . 

In the vertex model jTl, a rigidity transition takes place at 
a critical value of po = Po ~ When po < p^, cortical 
tension dominates over cell-cell adhesion and the energy bar¬ 
riers for local cell rearrangement and motion are hnite. The 
tissue then behaves as a elastic solid with hnite shear mod¬ 
ulus. When po > p^, cell-cell adhesion dominates and the 
energy barriers for local rearrangements vanish, resulting in 
zero rigidity and fluid-like behavior. While the energy func¬ 
tional for cell-cell interactions is identical in the vertex and 
SPV models, the two are truly distinct: the local minimum en¬ 
ergy states of the vertex model are not guaranteed to be similar 
to a Voronoi tessellation of cell centers, although we do find 
them to be very similar in practice. Therefore, we are also 
interested in whether a rigidity transition in the SPV model 
coincides with the rigidity transition of the vertex model. 

We define the effective mechanical interaction force ex¬ 
perienced by cell i as F, = —ViE (see Appendix [a| for de¬ 
tails). In contrast to particle-based models, F, is non-local 
and non-additive; F,- cannot be expressed as a sum of pairwise 
force between cells i and its neighboring cells. Nevertheless, 
one can show that momentum is still precisely conserved by 






























3 


this energy functional in the absence of the additional self¬ 
propulsion forces introduced below. 

In addition to F,, cells can also move due to self-propelled 
motility. Just as in SPP models, we assign a polarity vector 
hi = (cos 9,-, sin0,) to each cell; along hi the cell exerts a self¬ 
propulsion force with constant magnitude vo/p, where p is 
the mobility (the inverse of a frictional drag). Together these 
forces control the over-damped equation of motion of the cell 
centers r, 

dri 

— =pF, + voni. (2) 

dt 

The polarity is a minimal representation of the front/rear 
characterization of a motile cell Il50ll . While the precise mech¬ 
anism for polarization in cell motility is an area of intense 
study, here we model its dynamics as a unit vector that under¬ 
goes random rotational diffusion 

a,e, = Tii(0 

where 0, is the polarity angle that defines n,, and r|,(f) is a 
white noise process with zero mean and variance 2Dr. The 
value of angular noise Dr determines the memory of stochastic 
noise in the system, giving rise to a persistence time scale T = 
I/Dr for the polarization vector n. For small <C 1, the 
dynamics of n is more persistent than the dynamics of the 
cell position. At large values of Dr, i.e. when I/Dr becomes 
the shortest timescale in the model, Eq. (|^ approaches simple 
Brownian motion. 

The model can be non-dimensionalized by expressing all 
lengths in units of y/Ao and time in units of 1 /{iiKaAq). There 
are three remaining independent model parameters; the self 
propulsion speed vq, the cell shape index pQ, and the rota¬ 
tional noise strength Dr- We simulate a confluent tissue un¬ 
der periodic boundary conditions with constant number of 
N = 400 cells (no cell divisions or apoptosis) and assume that 
the average cell area coincides with the prefetTed cell area, 
i.e. (Aj) = Aq. This approximates a large confluent tissue in 
the absence of strong confinement. We numerically simulate 
the model using molecular dynamics by performing 10^ inte¬ 
gration steps at step size At = 10^^ using Euler’s method. A 
detailed description of the SPV implementation can be found 
in the Appendix Sec. 

II. CHARACTERIZING GLASSY BEHAVIOR 

We first characterize the dynamics of cell motion within the 
tissue by analyzing the mean-squared displacement (MSD) of 
cell trajectories. In Eig. [TJa), we plot the MSD as function 
of time, for tissues at various values of pQ and fixed vo = 0.1 
and Dr = 1. The MSD exhibits ballistic motion ( slope close 
to 2 on a log-log plot) at short times, and plateaus at inter¬ 
mediate timescales. The plateau is an indication that cells are 
becoming caged by their neighbors. Eor large values of po, 
the MSD eventually becomes diffusive (slope =1), but as po is 
decreased, the plateau persists for increasingly longer times. 


This indicates dynamical arrest due to caging effects and bro¬ 
ken ergodicity, which is a characteristic signature of glassy 
dynamics. 

Another standard method for quantifying glassy dynamics 
is the self-intermediate scattering function ll5^ : Fs{k,t) = 

giTAr(f) ^ Glassy systems possess a broad range of relax¬ 
ation timescales, which show up as a long plateau in Fs{t) 
when it is analyzed at a lengthscale q comparable to the near¬ 
est neighbor distance. Eig [T] (b) illustrates precisely this be¬ 
havior in the SPV model, when \k\ = k/tq, where ro is the 
position of the first peak in the pair cori'elation function. The 
average (...) is taken temporally as well as over angles of k. 
Fs{t) also clearly indicates that there is a glass transition as a 
function of po: at high po values F* approaches zero at long 
times, indicating that the structure is changing and the tissue 
behaves as a viscoelastic liquid. At lower values of po, F^ re¬ 
mains large at all timescales, indicating that the strvicture is ar¬ 
rested and the tissue is a glassy solid. Pig[2(d) demonstrates 
that at the strvictural relaxation time, the cell displacements 
show collective behavior across large lengthscales suggesting 
strong dynamical heterogeneity. This is strongly reminiscent 
of the ‘swirl’ like collective motion seen in experiment in ep¬ 
ithelial monolayers ll^l3l ITSl 1^ l40l . 


A. A dynamical order parameter for the glass transition 

Although the phase space for this model is three dimen¬ 
sional, we now study the model at a fixed value of = 1. 

We then search for a dynamical order parameter that dis¬ 
tinguishes between the glassy and fluid states as a func¬ 
tion of the two remaining model parameters, (vojPo)- A 
candidate order parameter is the self-diffusivity D/. Dg = 
lim,_>oo(Ar(f)^)/(4f). Eor practicality, we calculate D^ using 
simulation runs of 10^ time steps, chosen to be much longer 
than the typical caging timescale in the fluid regime. We 
present the self-diffusivity in units of Do = Vq/{ 2Dr), which is 
the free diffusion constant of an isolated cell. D^ff = DjDo 
then serves as an accurate dynamical order parameter that dis¬ 
tinguishes a fluid state from a solid (glassy) state in the space 
of (vojPo)^ matching the regimes identified using the MSD 
and Fq. In Eig. the fluid region is characterized by a finite 
value of Dgff and Dgjf drops below a noise floor of ~ 10^^ 
as the glass transition is approached. In practice, we label 
materials with D^ff > 10^^ as fluids indicated by an orange 
dot, and those with Dgff < 10^^ as solids indicated by blue 
squares. Importantly, we And that the SPV model in the limit 
of zero cell motility shares a rigidity transition with the vertex 
model Q at po ~ 3.81, and that this rigidity transition con¬ 
trols a line of glass transitions at finite cell motilities. Typical 
cell tracks (Eig. clearly show caging behavior in the glassy 
solid phase. 
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FIG. 2. (A) Glassy phase diagram for confluent tissues as function of cell motility vo and target shape index po at fixed = 1. Blue data 
points correspond to solid-like tissue with vanishing D^ff, orange points correspond to flowing tissues (finite Deff). The dynamical glass 
transition boundary also coincides with the locations in phase space where the structural order parameter q = {p/^/a) = 3.81 (dashed line). 
In the solid phase, ^ ~ 3.81 and ^ > 3.81 in the fluid phase. (B) Instantaneous tissue snapshots show the difference in cell shape across the 
transition. Cell tracks also show dynamical arrest due to caging in the solid phase and diffusion in the fluid phase. 


B. Cell shape is a structural order parameter for the glass 
transition 

In glassy systems it can be difficult to experimentally dis¬ 
tinguish between a truly dynamically arrested state and a state 
with relaxation times longer than the experimental time win¬ 
dow. Similarly, in tissues it is experimentally challenging to 
quantify a glass transition through the measurement of a dy¬ 
namical quantity such as the diffusivity D*. Identifying a static 
quantity that directly probes the mechanical properties of a tis¬ 
sue would therefore be a powerful tool for experiments. Pu- 
liafito et al. have suggested that shape changes accompany 
dynamical arrest in proliferating tissues HD. Similarly, a 
structural signature based on cell shapes - the shape index 
q = {jpj^fd) - was previously shown to be an excellent or¬ 
der parameter for the confluent tissue rigidity transition in the 
vertex model lIMI . In a model where cells were not motile 
(vQ = 0) we found that when po < 3.813, q is constant ^ 3.81 
and when po > 3.81 q grows linearly with po. Quite sur¬ 
prisingly, we found that the prediction of q = 3.813 works 
perfectly in identifying a jamming transition in in-vitro ex¬ 
periments involving primary human tissues, where cells are 
clearly motile (vo 7 ^ 0) Il3^ . At that time, we did not under¬ 
stand why the vo = 0 theory worked so well for these tissues. 

The prediction of a solid-liquid transition in the SPV 
model presented here provides an explanation for this obser¬ 
vation.We find that q (which can be easily calculated in ex¬ 
periments or simulations from a snapshot) can be used as a 
structural order parameter for the glass transition for all val¬ 
ues of Vo, not just at vq = 0. Specifically, the boundary defined 
by ^ = 3.813, shown by the blue dashed line in Fig.j^A) co¬ 


incides extremely well with the glass transition line obtained 
using the dynamical order parameter, shown by the round and 
square data points. The insets to Fig. also illustrate typical 
cell shapes; cells are isotropic on average in the solid phase 
and anisotropic in the fluid phase. This highlights the fact that 
q can be used as a structural order parameter for the glass tran¬ 
sition at all cell motilities, providing a powerful new tool for 
analyzing tissue mechanics. 


III. A THREE-DIMENSIONAL JAMMING PHASE 
DIAGRAM EOR TISSUES 

Having studied the glass transition as function of vq and 
Po at a large value of Dr, we next investigate the full three- 
dimensional phase diagram by characterizing the effect of Dr 
on tissue mechanics and structure. Dr controls the persistence 
time T = \/Dr and persistence length or Peclet number Pe ~ 
Vo/Dr of cell trajectories; smaller values of Dr correspond to 
more persistent motion. 

In Fig. I^A), we show several 2D slices of the three di¬ 
mensional jamming boundary. Solid lines illustrate the phase 
transition line identified by the structural order parameter 
q' = 3.813 as function of vq and po fo r a l arge range of Dr val¬ 
ues (from 10^^ to 10^). (In Appendix \B 2\ we demonstrate that 
the structural transition line q — 3.813 matches the dynamical 
transition line for all studied values of Dr.) In contrast to re¬ 
sults for particulate matter 0 , this figure illustrates that the 
glass transition lines meet at a single point (po = 3.81) in the 
limit of vanishing cell motility, regardless of persistence. 

Fig. [^B) shows an orthogonal set of slices of the jamming 
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FIG. 3. (A) The glass transition in vq — po phase space shifts as the persistence time changes. Lines represent the glass transition identified by 
the structural order parameter q = 3.81. The phase boundary collapse to a single point at pj = 3.81, regardless of Dr, in the limit vq —0. (B) 
The glass transition in po — Dr phase space shifts as a function of vq (from top to bottom: vq = 0.02,0.08,0.14,0.2,0.26) For large vq there is 
a crossover in the behavior at Dr ~ /jK^Aq = 1, as discussed in the main text. (C) The phase boundary between solid and fluid as function of 
motility vo> persistence I/Dr and po which is tuned by cell-cell adhesion can be organized into a schematic 3D phase diagram. Red lines on 
the surface correspond to iso-VQ contours and blue lines correspond to iso-D^ contours. 


diagram, illustrating how the phase boundary shifts as func¬ 
tion of po and Dr at various values of vq- This highlights the 
interesting result that a solid-like material at high value of Dr 
can be made to flow simply by lowering its value of Dr- The 
crossover in behavior at large vq occurs when the persistence 
time I/Dr is approximately equal to the viscous relaxation 
time 1/{pKaAo) = 1. 

These slices can be combined to generate a three- 
dimensional jamming phase diagram for confluent biological 
tissues, shown in Fig.[^C). This diagram provides a concrete, 
quantifiable prediction for how macroscopic tissue mechanics 
depends on single-cell properties such as motile force, persis¬ 
tence, and the interfacial tension generated by adhesion and 
cortical tension. 

We note that Fig. [^C) is significantly different from the 
jamming phase diagram conjectured by Sadati et al ll42l . 
which was informed by results from adhesive particulate mat¬ 
ter ll55l . For example, in particulate matter adhesion enhances 
solidification, while in confluent models adhesion increases 
cell perimeters/surface area and enhances fluidization. In ad¬ 
dition, we identify “persistence” as a new axis with a poten¬ 
tially significant impact on cell migration rates in dense tis¬ 
sues. 

To better understand why persistence is so important in 
dense tissues, we first have to characterize the transitions be¬ 
tween different cellular structures. In the limit of zero cell 
motility, the system can be described by a potential energy 
landscape where each allowable arrangement of cell neigh¬ 
bors corresponds to a metastable minimum in the landscape. 
There are many possible pathways out of each metastable 
state: some of correspond to localized cell rearrangements, 
while others correspond to large-scale collective modes. The 
maximum energy required to transition out of a metastable 
state along each pathway is called an energy barrier 10. 

We observe that tissue fluidity can increase drastically with 
increasing Dr at finite cell speeds. This suggests that different 



FIG. 4. (A-C) Instantaneous cell displacements at po = 3.65 

and Vq = 0.5. They are different from the displacements shown in 
Fig .[TJD) which are averaged over the structural relaxation timescale. 
(A) At the lowest value of Dr = 0.01, the cells are able to flow by 
collectively displacing along the ‘soft’ modes of the system (Ap¬ 
pendix. |B l[ l. (B) At Dr = 0.1, collective displacements are less 
pronounced. (C) For Dr = I and larger, the displacements appear 
disordered and uncorrelated. 


pathways (with lower energy barriers) must become dynami¬ 
cally accessible at higher values of Dr- 

One hint about these pathways comes from the instanta¬ 
neous cell displacements, shown for different values of Dr in 
Fig.0 At high values of Dr, (po = 3.78, vq = 0.1) the instan¬ 
taneous displacement field is essentially random and largely 
uncorrelated, as shown in Fig. and the material is solid¬ 
like. There is no collective behavior among cells, and each 
cell ‘rattles’ independently near its equilibrium position. 

However, as Dr is lowered, the instantaneous displacement 
field becomes much more collective (Fig. and the tissue 
begins to flow, presumably because these collective displace¬ 
ment fields correspond to pathways with lower transition en¬ 
ergies. 

Two obvious questions remain: How does a lower value 
of Dr generate more collective instantaneous displacements? 
Why should collective instantaneous displacements generi- 
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cally have lower energy barriers? The first question can be 
answered by extending ideas first proposed by Henkes, Fily 
and Marchetti 1201 to explain why motion in self-propelled 
particle models seems to follow the ‘soft modes’ of a solid. 
This argument is based on a simple, yet powerful observation: 
in the limit of zero motility (vo = 0), a solid-like state will have 
a well-defined set of normal modes of vibration (with frequen¬ 
cies {cOv}), and a corresponding set of eigenvectors ({cv}) that 
forms a complete basis. At higher motilities (vq > 0) near the 
glass transition, the motion of particles in the system can be 
expanded in terms of the eigenvectors. As discussed in Ap¬ 
pendix |B^ one can use this observation to show that in the 
limit of Dr 0, motion along the lowest frequency eigen- 
modes is amplified - the amplitude along each mode is pro¬ 
portional to 1/(0^. These low-frequency normal modes are 
precisely the collective displacements observed for low D, . 

The second question is more difficult to answer because it 
is impossible to enumerate all of the possible transition path¬ 
ways and energy barriers in a disordered material. However, 
a partial answer comes from recent work in disordered par¬ 
ticulate matter showing that low-frequency normal modes do 
have significantly lower energy barriers 1^ IfiTl than higher 
frequency normal modes. A similar analysis could potentially 
be performed in vertex or SPV models. 

IV. A CONTINUUM MODEL FOR GLASS TRANSITIONS 
IN TISSUES 

Although continuum hydrodynamic equations of motion 
have been developed by coarse graining SPP models in the 
dilute limit, there is no existing continuum model for a dense 
active matter system near a glass transition. Here we propose 
that a simple trap or Soft Glassy Rheology (SGR) 1471 
model provides an excellent continuum approximation for the 
phase behavior in the large Dr Brownian regime, but fails in 
the small Dr limit. 

For large Dr it is known that particle behave like Brownian 
particles with an effective temperature Teff = vl/2^lDr mi. 
This mapping becomes exact when at fixed “effec¬ 

tive inertia” {^iDr)^^ lfT3l . In other words, like in granular 
systems ffua, the effective temperature in SPP is dominated 
by kinetic effects. Guided by this result we conjecture that 
in our model the temperature also scale quadratically with the 
velocity, 

TeZ/occVo. (4) 

Physically, this effective temperature gives the amount of en¬ 
ergy available for individual cells to vibrate within their cage 
or ‘trap’. 

The next important question is how to characterize the ‘trap 
depths’, or energy barriers between metastable states. In the 
Brownian regime (large Dr) there is no dynamical mechanism 
for the cells to organize collectively, and therefore a reason¬ 
able assumption is that the rearrangements are small and lo¬ 
calized. 

In Is), some of us explicitly calculated the statistics of en¬ 
ergy barriers for localized rearrangements in the equilibrium 



FIG. 5. Comparison between SPV glass transition and an analytic 
prediction based on a Soft Glass Rheology (SGR) continuum model. 
The dashed line corresponds to an SGR prediction with no fit param¬ 
eter based on previously measured vertex model trap depths O. Data 
points correspond to SPV simulations with Dr = 10^^ and where we 
have defined T^ff = cvg with c = 0.1 as the best-fit normalization pa¬ 
rameter. Blue points correspond to simulations which are solid-like, 
with Dgff < 10^^, and the boundary of these points define the ob¬ 
served SPV glass transition line. (Inset) L 2 difference between SPV 
glass transition line (at best-fit value of c) and the predicted SGR 
transition line at various values of Dr- The SGR prediction based on 
localized T1 trap depths works well in the high Dr limit, but not in 
the low Dr limit. 


vertex model. In the 2D vertex model, one can show that 
localized rearrangements must occur via so-called T1 tran¬ 
sitions ll58l . Using a trap model lf3^ or Soft Glassy Rheol- 
ogy(SGR) BtII framework, we were able to use these statis¬ 
tics to generate an analytic prediction, with no fit parameters, 
for the glass transition temperature Tg as function of po- 

To see if the SGR prediction for the glass transition holds 
for the SPV model in the large Dr limit, we simply overlay 
the data points corresponding to glassy states from the SPV 
model with the glass transition Tg line predicted in IH. There 
is one fitting parameter c that characterizes the proportion¬ 
ality constant in Eq. Fig shows that the SPV data for 
Dr = 10^ is in excellent agreement with our previous SGR 
prediction. Because Teff ~ vl, and the glass transition line 
scales as Tg ^ p^— po for po p^ and Dr -> 0, the glass 
transition line scales as vq ~ (pg — po) in those limits. 

The reason the effective temperature SGR model works 
here is that, like in SPP models of spherical active Brownian 
colloids, the angular dynamics of each cell evolves indepen¬ 
dently of cell-cell interactions and of the angular dynamics of 
other cells. An additional alignment interaction that couples 
the angular and translational dynamics may therefore modify 
this behavior. 

To our knowledge, this is the first time that a SGR / trap 
model prediction has been precisely tested in any glassy sys¬ 
tem. This is because, unlike most glass models, we can enu¬ 
merate all of the trap depths for localized transition paths in 
the vertex model. 
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However, for small values of Dr, we have shown that cell 
displacements are dominated by collective normal modes, and 
therefore the energy barriers for localized T1 transitions are 
probably irrelevant in this regime. The inset to Fig shows 
the deviation (L^-norm) between glass transition lines in the 
SPV model and Tl-based SGR prediction as a function of Dr. 
We see that the SGR prediction fails in the small Dr limit, 
as expected. A better understanding of the energy barriers 
associated with collective modes will be required to modify 
the theory at small Dr- 

V. DISCUSSION AND CONCLUSIONS 

We have shown that a minimal model for confluent tissues 
with cell motility exhibits glassy dynamics at constant den¬ 
sity. This model allows us to make a quantitative prediction 
for how the fluid-to-solid/jamming transition in biological tis¬ 
sues depends on parameters such as the cell motile speed, the 
persistence time associated with directed cell motion, and the 
mechanical properties of the cell (governed by adhesion and 
cortical tension). We define a simple, experimentally acces¬ 
sible structural order parameter - the cell shape index - that 
robustly identifies the jamming transition, and we show that a 
simple analytic model based on localized T1 rearrangements 
precisely predicts the jamming transition in the large Dr limit. 
We also show that this prediction fails in the small Dr limit, 
because the instantaneous particle displacements are domi¬ 
nated by collective normal modes. 

This model makes several experimentally-verifiable predic¬ 
tions for cell shape and tissue mechanics: 

• The order parameter q' = 3.81 is a structural signature 
for the glass transition, even in tissues with significant 
cell motility or dynamics. This prediction has already 
been tested in epithelial lung tissue lf38l . but it should 
be much more broadly applicable. We have performed 
a rudimentary shape analysis of a small number of im¬ 
ages from other systems that have been previously pub¬ 
lished, including proliferating MCDK monolayers cni 
and convergent extension in fruit fly development im 
and found that the shapes are consistent with this pre¬ 
diction. A much more careful analysis with full data 
sets should be performed to further validate this predic¬ 
tion or understand where it breaks down. 

• In the limit of vanishing cell motility, shape and pres¬ 
sure fluctuations should vanish when the jamming tran¬ 
sition is approached from the solid side, and remain 
zero in the fluid. A finite motility vq will induce such 
fluctuations in the fluid phase, as confirmed by prelimi¬ 
nary calculation of cellular stresses and pressure in the 
SPV model ll62ll . This could be studied by combin¬ 
ing measurement of cell shape fluctuations with traction 
force microscopy (TFM) in wound healing assays. Af¬ 
ter locating the glass transition by imaging cell shape 
changes, it may be possible to extract information on 
cell motility vq from cellular stresses and pressure in¬ 
ferred from TFM in the fluid phase near the glass tran¬ 


sition. This suggests that one may estimate cell motility 
by examining the changes in cellular stresses and pres¬ 
sure in the cell monolayer near the unjamming transi¬ 
tion and assuming that the local velocity of the mono- 
layer is very small just above the transition. The latter 
assumption can also be verified independently via par¬ 
ticle image velocimetry (PIV). 

• Cell proliferation, so far neglected in our model, causes 
an increase in cell number density in confluent tissues. 
Often this is accompanied by a reduction in individual 
cell motility vo, via contact inhibition of locomotion. In 
cases where this is the dominant effect and changes to 
the ratio between Aq and Pq negligible, our work 
predicts that proliferation would drive the system to¬ 
wards jamming. This is consistent with existing reports 
in the literature ifTOll . although more work is required to 
test the prediction carefully. In tissues where vq remains 
low at all times ED, our model predicts that prolifera¬ 
tion can either cause jamming or unjamming, depend¬ 
ing on whether cell divisions are oriented in such a way 
to decrease or increase cell shape anisotropy. 

• Spatial correlations and fluctuations of the cell dis¬ 
placement field, such as swirl sizes (0012312911401? 
, should grow as a tissue approaches the glass transi¬ 
tion from the fluid side. Very recently, a similar pre¬ 
diction for displacements and correlation lengths based 
on a particle-based model has been verified in one cell 
type im. The SPV model, which makes predictions for 
cell shapes in addition to displacements and correlation 
lengths, could be tested simply by compiling detailed 
statistics about cell shapes and cell motion in epithelial 
monolayers. 

Although all of the work presented here focuses on the SPV 
model, which tracks cell centers and therefore has only two 
degrees of freedom per cell, we found that in the limit of zero 
cell motility it exhibits the same rigidity transition as the ver¬ 
tex model which has two degrees of freedom per vertex. We 
have also checked that an ’’active vertex model”, where ac¬ 
tive motile forces are added to the vertex model vertices, also 
exhibits a robust glass transition characterized by the shape 
order parameter q. The fact that two models with ostensibly 
different degrees of freedom share the same transition sug¬ 
gests that there is a deeper universality, perhaps generated by 
isostaticity, that remains to be understood. 

Another result of this work is the surprising and unexpected 
differences between confluent models (such as the vertex and 
SPV models) and particle-based models (such as Lennard- 
Jones glasses and SPP models). For example, works by 
Berthier ||6l, Fily and Marchetti lfT4l in SPP models suggest 
that the location of the zero motility glass transition packing 
density (|)c (defined as the density at which dynamics cease 
in the limit of vq —>^ 0) depends the value of noise. Dr. This 
is also related to the observation that the jamming and glass 
transition are not controlled by the same critical point in non¬ 
active systems EH [371. We And this is not the case in the SPV 
model. Fig.j^a & b) show that while the glass transition point 
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Pq shifts with Dr at finite values of vq, in the limit of vanishing 
motility, all glass transition lines merge on to a single point in 
the limit vq —>■ 0, namely /7 q = 3.81. 

Given these differences, it is important to ask which type 
of model is appropriate for a given system. We argue that 
SPV models are maybe more appropriate for many biologi¬ 
cal tissues. Whereas SPP models interact with two-body in¬ 
teractions that only depend on particle center positions, both 
SPV and vertex models naturally incorporate contractility as 
a key property of living cells and capture the inherently multi¬ 
body nature of intercellular forces due to shape deformations. 
Unlike equilibrium vertex models, SPV models account for 
cell motility, and they are also much easier to simulate in 3D 
(which is nearly impossible in practice for the vertex model.) 

Recent work by Li and Sun 12^ also models a confluent 
cell as a Voronoi tessellation of the plane. An important dif¬ 
ference between this work and ours is that in Ref. ESil cell¬ 
cell adhesion is captured via a potential that is quadratic in the 
distance between cell centers, just as in particle models. We 
might guess that stronger cell-cell adhesion in their model will 
result in stiffening of the tissue, which is common for parti¬ 
cle based models, although that remains to be tested in active 
systems. In contrast, adhesion enters our model through the 
coupling of the shape energy to the cell perimeter. Increas¬ 
ing cell-cell adhesion (or decreasing cortical tension) yields a 
larger value of pQ, which leads to the tissue becoming softer. 

We expect that other shape-based models of confluent tis¬ 
sue dynamics will also yield the glass transition described 
here. For example, it has been reported in recent works based 
on the Cellular Potts model Eiiia and in a modified SPP 
model ca when the cell motile force is decreased beyond a 
certain threshold, the motion of cells transitions from diffusive 
to sub-diffusive. This is similar to crossing the glass transition 
line in the SPV model by decreasing the value of vq. 

In this work and in previous work based on the vertex 
model, the cell volume is generally assumed to be fixed. 
While this is a good assumption in developmental systems 
such as drosopholia Qaiia and zebrahsh ED, epithelial 
tissues cells can show significant volume fluctuations, as re¬ 
ported recently MM- Therefore, it will be important to 
incorporate volume fluctuations in future iterations of the ver¬ 
tex model or the SPV model, as they introduce another source 
of active shape fluctuation and could therefore lead to jam¬ 
ming or unjamming of the tissue locally and potentially shift 
the location of the rigidity and glass transitions. 

In our version of the SPV model, we have assumed that cell 
polarity is controlled by simple rotational white noise. It is 
also possible to include more complex mechanisms. For ex¬ 
ample, external chemical or mechanical cues could be mod¬ 
eled by coupling vq and n, to chemoattractant or mechanical 
gradients, allowing waves or other pattern formation mecha¬ 
nisms to interact with the jamming transition. Similarly, sim¬ 
ple alignment rules (such as those in the Viscek model Ell) 
could lead to collective flocking modes that also affect glassy 
dynamics. 

Another interesting extension of the SPV model would be 
to study the role of cell-cell friction, which has already been 
shown to be important in controlling collective dynamics in 


particle-based tissue models Ea. Our current model includes 
viscous frictional coupling of cell to the 2D substrate and cell¬ 
cell adhesion enters as a negative line tension on interfaces. 
However, it would be possible to add a frictional force be¬ 
tween cells proportional to the length of the edge shared be¬ 
tween two cells, and we know from previous work on par¬ 
ticulate glasses that these localized frictions can change the 
location of jamming/glass transition and the nature of spatial 
correlations in a glass im sa¬ 
lt is also tempting to speculate about the relationship be¬ 
tween the unjamming transition captured by our model and the 
epithelial-mesenchimal transition (EMT) that precedes cell 
escape from a solid tumor mass. The EMT involves signifi¬ 
cant changes in cell-cell adhesion and cytoskeletal composi¬ 
tion, with associated changes in cell shape and motility. This 
suggests that escape from the tumor mass is controlled not 
just by the chemical breakdown of the basement membrane, 
but also by specific changes in mechanical properties of both 
individual cells and the surrounding tissue ll60ll . One could 
then hypothesize that the collective unjamming described here 
may provide the first necessary step towards the mechanical 
changes needed for cell escape from primary tumors. 

In particular, recent work suggests that cancer tumors are 
mechanically heterogeneous, with mixtures of stiff and soft 
cells that have varying degrees of active contractility E9l . Our 
jamming phase diagram suggests that the soft cells, which of¬ 
ten exhibit mesenchymal markers and presumably correspond 
to higher values of po, might unjam and move towards the 
boundary of a primary tumor more easily than their stiff coun¬ 
terparts. Examining the effects of tissue heterogeneity on tis¬ 
sue rigidity and patterns of cell motility is therefore a very 
promising avenue for developing predictive theories for tumor 
invasiveness and metastasis. 


Appendix A: Simulation algorithm for the SPV model 

To create an initial configuration for the simulation, we hrst 
generate a seed point pattern using random sequential addi¬ 
tion (RS A) and anneal it by integrating Eq. |^with vq = 0 
for 100 MD steps. The resulting structure then serves as an 
initially state for all simulations runs. The use of (RSA) only 
serves to speed up the initial seed generation as using a Pois¬ 
son random point pattern does not change the results presented 
in this paper. 

At each time step of the simulation, a Voronoi tesselation is 
created based on the cell centers. The intercellular forces are 
then calculated based on shapes and topologies of the Voronoi 
cells (see discussion below). We employ Euler’s method to 
carry out the numerical integration of Eq.|^ i.e., at each time 
step of the simulation the intercellular forces is calculated 
based on the cell center positions in the previous time step. 

In a Delaunay triangulation, a trio of neighboring Voronoi 
centers define a vertex of a Voronoi polygon. Eor example in 
Fig. s (r,,rj,f>) define the vertex h$, which is given by 

^3 = an + Pg + y4, (A1) 
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where the coefficients are given by 

a = \\rj - hfin - rj) ■ in - rk)/D 
P = \\h - rkfin - ^i) ■ in - rk)/D 
y=\\n-rjf{h-h)-{h-n)/D 
D = 2\\{ri-rj) x {rj-fk)\\^. 



FIG. 6. Cell centers positions are specified by vectors {?}. They 
form a Delaunay triangulation (black lines). Its dual is the Voronoi 
tessellation (red lines), with vertices given by {/i}. 


In the vertex model, the total mechanical energy of a tissue 
depends only on the areas and perimeters of cells: 

N 

£ = ^ [KpiA^-Ao)^+Kp{Pt-Po)^] . (A3) 

!=1 

In a Voronoi tessellation, the area and perimeter of a cell i can 
be calculated in terms of the vertex positions 


since the interaction is inherently multi-cellular in nature and 
interactions between i and / also depend on k and / (see 

Fig-0- 

For the typical configuration shown in Fig.0 the first term 
in Eq.[^ can be expanded using the chain rule and calculated 
using Eq. 


dE^ 

drip 


L 


/ dEj dh2v 
\dh2v drt^ 


dEj dhjiM 
dh^v drifj 


(A7) 


In Eq. |A7[ only terms involving h 2 and are kept since Ej 
does not depend on other vertices of cell i. v is a cartesian 
coordinate index. The energy derivative in Eq. A7 can be cal¬ 


culated in a straightforward way, by using Eqs. A3 and A4 


dE^ 

dh2x 


=2Ka(Aj - Ao) ^ + 2Kp{Pj - Fo) ^ 
dh2x dh2x 

=Ka {A-j — AQ){h2y — h^y) 


P2Kp{Pj-PQ) 


^2x ^Ix ^ ^2x hjfx I 

11^7-^211 11^2-^311/ 


(A8) 


and 


dE^ 

dh2y 


=2KAAj-A«)^+2KAP,-P„)^ 

=Ka (a j—Ao) {hix — hx) 


+ 2Kp{Pj-Po) 


( ^2v ^7y 

Vll^7-^2|| 


^2y h-yy I 

\\h2-%\\j 


(A9) 


Similarly, the second term on the r.h.s. of Eq. [A5| can be cal¬ 
culated in a similar way. 


7i-l 

Pi — ll^m 

m=0 

— 2 11 ^ ^m+1 11 j 

^ m=0 


Appendix B: Cell displacements and structnral order parameter 
as a fnnction of Dr 

(A4) 

1. Expanding cell displacements in an eigenbasis associated 
with the nnderlying dynamical matrix 


where Zi is the number of vertices for cell i (also number of 
neighboring cells) and m indexes the vertices. We use the con¬ 
vention h^. = ho. 

With these definitions, the total force on cell-/ can be cal¬ 


culated using Eq. A3 


- 

dn,, 


3k- 

j€n.n.{i) 


dEi 

3r 


(A5) 


here p denotes the cartesian coordinates (x,y). The first term 
on the r.h.s. of Eq. A5 sums over all nearest neighbors of cell 
i. It is the force on cell / due to changes in neighboring cell 
shapes. The second term is the force on cell / due to shape 
changes brought on by its own motion. 

dE 

It maybe tempting to treat ^ as the force between cells-/ 
and j, but 


dEj ^ dEi 
drip drjy, 


(A6) 


In the absence of activity (vq = 0), the tissue is a solid 
for po < Pq = 3.81. As Vo is increased, the solid behavior 
persists up to vq = Vq(po), which is given by the glass tran¬ 
sition line in Eig. 2. In order for the tissue to flow, suffi¬ 
cient energy input is needed to overcome energy barriers in 
the potential energy landscape, which are a property of the 
underlaying solid state at vo = 0. In this limit, the instan¬ 
taneous cell center positions {ri{t)} can be thought of as a 
small displacement {di{t)} from the nearest solid reference 
state {fpi} EOl where di{t) = Vi — ro,. The /p,- correspond 
to positions of cell in a solid, which has a well-defined lin¬ 
ear response regime Q. The linear response is most con¬ 
veniently expressed as the eigen-spectrum of the dynamical 
matrix T>,;a|3- Since the eigenvectors {e,',v} of Dyap form a 
complete orthonormal basis, the cell center displacement can 
then be expressed as a linear combination of {e,,v} 

di{t) =Y^ay{t)eiy (Bl) 

V 
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For simplicity, we will adopt the Bra-ket notation and express 
the eigenbasis simply as |v) and Eg. |B1 [ becomes 

M>=L«v(0|v), (B2) 

V 

where 

O|v) = (o2|v) (B3) 

and cOy are the eigenvalues of the dynamical matrix. 

The polarization vector n, can also be expressed as a linear 
combination of eigenvectors 

|n)=^^(f)|v). (B4) 

V 

Since the polarization vector and eigenvector are both unit 
vectors, it follows that by,{t) = (n|v) = coi(0v — \(/). Where \(/ 
is the angle of the polarization and 0v the angle of the eigen¬ 
vector. 

Then the equation of motion for di (Eg. 2), can be rewritten 
as 


to get the ensemble averaged solution for the amplitude be¬ 
comes 

(av(0) = av(0)e^^' + vocos (0v - \(/(0)) ^—■. 

Dr — k 

(BID 

In the limit of Dr ^ and Eg. |B1 1 [ becomes 

(flv(O) =«v(0)e^'^'. (B12) 


This suggests that while normal modes control the rate of de¬ 
cay, they do no affect the long-time behavior. 

However as Dr —^ 0, Eg. B1 l[becomes 


av(f) = ^ cos (0v - ^(0)) (l - 

(B13) 

The second term in this equation scales as ^ 1 /cOy. Therefore, 
at short times (corresponding to instantaneous response), the 
mode amplitude ay is much larger for modes at lower frequen¬ 
cies. Since the reference state is an elastic solid with Debye 
scaling D{(o) ^ co as (O —0 ||71, this suggests that the dis¬ 
placement will be heavily dominated by the lowest frequency 
modes that are spatially more collective in nature. 


d = —p 


dE 

dri 


- von,' 


(B5) 


Using Egs. [^ - B5 we find 

J^{v\d) = -^l{v\b d)+VQ{y\n) ,or 
^av(D = -pCOvav(f)+ ''’0^v(f)- 

Then the equation of motion for each amplitude is 
d 


(B6) 


dt 


ay{t) = —pa)yav(D + ''ocos(0v —1|/) 

¥ = T1- 


(B7) 


2. Effect of Dr on glass transition boundary 


Dr = 0.01 


Dr = too 




1.0 

0.8 

0.6 

D.4 

0.2 

0.0 


3.0 3.2 3.4 3.6 3.8 4.0 
PO 


FIG. 7. Comparison between glass transition boundaries obtained 
using shape order parameter (red line) and £>e// (blue squares and 
orange circles). 


This is just the equation of motion for a self-propelled particle 
tethered to a spring with active forcing that is strongest along 
the direction of the eigenvector 0- The solution is then: 

ay(t) = ay(t = + V0 f dt'e^^^‘^’'^cos(Qy — \\t), (B8) 

JO 

where k = pco^. 

Solving for the ensemble averaged quantity: 

(av(D) = «v(f = + Vo [ ^(coi(0v — 

JO 

(B9) 

and using the relations 

(cost|/(f)) = cos\(/(0)e^^'''; 

(sint|/(f)) = sint(/(0)e^^''' (BIO) 

cos(0v —'ll) = sin(0v)sin(\(/) -|-cos(0v)cos(t)/), 


Figure. [^ shows that the location in phase space where the 
shape index q' = 3.81 is in excellent agreement with the dy¬ 
namical solid-fluid phase boundary determined by Defj, at all 
values of Dr- 
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